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We present a model for the growth of the transportation network inside nests of the social insect 
subfamily Termitinae (Isoptera, termitidae). These nests consist of large chambers (nodes) con¬ 
nected by tunnels (edges). The model based on the empirical analysis of the real nest networks 
combined with pruning (edge removal, either random or weighted by betweenness centrality) and 
a memory effect (preferential growth from the latest added chambers) successfully predicts emer¬ 
gent nest properties (degree distribution, size of the largest connected component, average path 
lengths, backbone link ratios, and local graph redundancy). The two pruning alternatives can be 
associated with different genuses in the subfamily. A sensitivity analysis on the pruning and mem¬ 
ory parameters indicates that Termitinae networks favor fast internal transportation over efficient 
defense strategies against ant predators. Our results provide an example of how complex network 
organization and efficient network properties can be generated from simple building rules based on 
local interactions and contribute to our understanding of the mechanisms that come into play for 
the formation of termite networks and of biological transportation networks in general. 

PACS numbers: 89.75.He, 89.75.Fb, 05.65.+b, 87.23.Cc 


I. INTRODUCTION 

Biological transportation networks are a fundamental 
component of all living systems, allowing the exchange 
of information and material at the scale of the whole sys¬ 
tem. It is unsurprising that structures specialized for 
transportation are found at all levels of biological orga¬ 
nization, including intracellular transportation along the 
cytoskeleton, and vascular and neural networks within 
the body of animals. At the biological scale of ani¬ 
mal groups and communities, the transportation of food 
and material, the movements of animals and the inter¬ 
individual encounters are often supported by specialized 
transportation networks of trails, galleries, and burrows 
that the animals produce and use in their exploration and 
foraging movements (reviewed in 0 )- Social insects in 
particular are known to produce some of the most com¬ 
plex networks of trails and galleries in the animal king¬ 
dom, probably as a result of their high level of sociality. 
These networks include trail networks in ants or termites 
2r4J, and underground systems of tunnels formed by ants 
3 or termites [a-Q. 

All biological transportation networks share similar 
functions: by favoring transportation over distances 
much larger than those permitted by simple diffusion 
they mediate the integration of the different parts that 
compose a biological system, supporting the functional 


unity of the system as a whole. They also share a sim¬ 
ilar morphogenesis, in the sense that almost all biolog¬ 
ical transportation networks are produced as the result 
of self-organized (SO) morphogenetic processes whereby 
the growth is driven by locally available information, in 
the absence of a pre-existing master plan of the net¬ 
work 00- The formation of animal and human trans¬ 
portation networks has been modeled with models based 
on growth alone (e. g., a rmy ant raid networks EH , ant 
gallery networks 0 [l2[| , or urban street patterns jl3l|). 
pruning of an existing network , or a combination 

of both [Ml- 

A common characteristic of all these networks is that 
they are embedded in a 2D or 3D environment, that is, 
both the position of network nodes and the layout of net¬ 
work edges are associated with sets of spatial coordinates. 
The effects of spatial embedding cannot be neglected 
when trying to understand the formation and topological 
properties of these networks, because the probability of 
existence of a connection between two nodes depends not 
only on their relative distance [19| , but also on the phys¬ 
ical arrangement and steric interactions between edges. 
While a large part of the existing network literature has 
dealt with social or communication networks, which are 
comparatively less affected by spatial constraints (see re¬ 
views in [20l - l22f ). the theoretical foundations underlying 
the analysis and the modeling of spatially explicit net- 
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works are comparatively less developed (see 0 for a 
review) and often deal with specific fields such as urban 
transportation and human mobility patterns [24| . 

In this paper we focus on a specific class of biological 
transportation networks represented by the network of 
chambers and tunnels that termites of the subfamily Ter- 
mitinae produce in the above-ground part of their nest. 
We identify how structural features of these transporta¬ 
tion networks emerge through self-organization based on 
local rules. The nodes of these networks correspond 
to chambers and the edges to tunnels connecting these 
chambers. These networks provide safe living space to 
termites and are connected to an underground tunnel 
network (not considered in this study) that connect the 
nest to foraging grounds 0 - 127} . Termite movements 
along these networks involve both bringing back food and 
distributing it to the colony (e.g. termite soldiers, larvae 
and the royal couple) and daily patterns of movements 
of various individuals to different parts of the nest with 
favorable environmental conditions. 

By analyzing an extended dataset containing nests 
from three Termitinae genuses: Cubitermes , Procubiter¬ 
mes and Thoracotermes , we formulate a model for the 
growth of the transportation networks internal to the 
nests of these different termites. The model has two 
free parameters, A that controls preferential growth on 
the periphery and £ that determines the pruning or re¬ 
moval of existing edges. Two variants of pruning, random 
pruning or pruning weighted by betweenness centrality, 
are explored, and we compare both variants to random 
geometric graphs. We calibrate A and £ to the real nests 
and show that network properties not used in the cali¬ 
bration process are faithfully reproduced by the models 
with pruning, the variant based on betweenness central¬ 
ity performing slightly better. The validated models are 
then used to assess the sensitivity of some network prop¬ 
erties, linked to internal transportation or nest defense, 
to the two free parameters in order to interpret the values 
that have been found for the real nests. 


II. MATERIALS AND METHODS: FROM AN 
EMPIRICAL NEST ANALYSIS TO THE NEST 
GROWTH MODEL 

A. The available termite nests 

The 12 networks modeled in this paper (see Table |]J 
correspond to the above ground part of termite nests 
from the African continent. They all belong to the Ter¬ 
mitinae sub-family that have a common architecture: dis¬ 
tinct chambers with a large diameter interconnected by 
tunnels of a small diameter (these different elements can 
be easily identified because in these nests the connect¬ 
ing tunnels have a diameter about 10 times smaller than 
the chambers’ diameter). In our network representation, 
chamber barycentres are associated to nodes and tunnels 
between chambers to edges between the corresponding 


nodes. The detailed network extraction method from x- 
ray tomographies of the nests has been described in 0 . 
Node coordinates are converted from voxel positions to 
metric (x, y , z) coordinates. The z-axis corresponds to 
the vertical direction. 

The majority of the nests used in our study have been 
obtained from natural history museums in France, with 
their taxonomic identity only known to the genus level. 
There are six Cubitermes nest networks (already pub¬ 
lished in [28}), four Procubitermes sjoestedti networks col¬ 
lected in 2007 in Cote d’Ivoire, and two Thoracotermes 
networks (the larger one, a T. macrothorax , was collected 
in the Republic of the Congo in 2009). See the SM 0 
for the general shape of all nests. While some aspects 
of mound architecture are typical of nests of each genus 
(mushroom like shapes in Cubitermes , straight pillars in 
Thoracotermes , and conic forms in Procubitermes ), at¬ 
tempts to use this architecture for taxonomy have failed 
0 ] , pointing to the fact that all the nests in the Termiti¬ 
nae subfamily (hence all the nests analyzed here) share 
similar morphological characteristics and possibly result 
from similar construction rules. 


B. Empirical network analysis 

Each tunnel or edge in the nest networks can be con¬ 
sidered as a three dimensional vector represented by 
Cartesian coordinates ( x , y , z) or by spherical coordinates 
(r, 0, </>). The latter will be used to characterize the nests. 
Since there is no natural orientation in a tunnel connect¬ 
ing chambers = (xi,yi,z\) and 17 ’2 = (% 2 ,y 2 , Z 2 ), it 
can be represented as either c \ — 2 or as c 2 — c 1 

(Fig. OJa)). Both vectors are used in the empirical de¬ 
scription of the nests, e.g., FigOJb-d) for nes t C9 that 
shows the distributions of tunnel length r, the vertical 
component 0 and the horizontal component (j>. The dis¬ 
tributions for the other 11 nests are shown in the SM 



Both the r distribution (Fig. [Ub)) and the 9 distribu¬ 
tion (Fig.QIc)) resemble normal distributions, the latter 
with mean ir/2. However, the distribution of <f> (x — y 
plane) rather follows a uniform distribution (Fig. [l](d)). 
These patterns are confirmed in the other nests (see SM 
0). This means that each termite nest can be charac¬ 
terized by three parameters: the mean and standard de¬ 
viation of its tunnel length distribution (f, ay), and the 
standard deviation of the vertical 9 direction component 

Oe)- 

C. Model description 

Based on the above empirical observations we propose 
a simple nest growth model, betweenness based pruning 
(BBP) model, with five main procedures: (i) set nest 
boundaries and initial node, (ii) determine the initiation 
node for the next tunnel, (iii) construct an edge (iv) ere- 
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TABLE I. General statistics of the analyzed termite nest networks and the model parameters that have been calibrated to each 
nest (see model description). Nr-C: number of chambers; Nr-T: number of tunnels; ND: average node degree. LCC: size of the 
largest connected component. (A bbp,£,bbp) are the parameters for the BBP model, (A rp,^rp) the ones for the RP model, 
and Rrgg is the parameter for the RGG model. 


Nest 

Genus 

Nr-C Nr-T 

ND 

LCC A bbp £bbp 

A RP 

£rp 

Rrgg 

C9 

Cubitermes 

532 

682 

2.56 

507 

0.021 

0.117 

0.020 

0.126 

0.199 

CIO 

Cubitermes 

396 

371 

1.87 

349 

0.065 

0.292 

0.065 

0.292 

0.197 

Cll 

Cubitermes 

344 

310 

1.80 

260 

0.066 

0.279 

0.065 

0.274 

0.205 

C12 

Cubitermes 

190 

234 

2.46 

183 

0.067 

0.246 

0.066 

0.244 

0.274 

C18 

Cubitermes 

312 

343 

2.2 

287 

0.031 

0.250 

0.029 

0.248 

0.226 

C19 

Cubitermes 

295 

445 

3.02 

268 

0.045 

0.085 

0.046 

0.080 

0.258 

P67 

Procubitermes 

1123 

2149 

3.83 

1091 

0.025 

0.0 

0.025 

0.0 

0.177 

P78 

Procubitermes 

675 

878 

2.60 

598 

0.014 

0.046 

0.014 

0.044 

0.185 

P79a 

Procubitermes 

440 

525 

2.39 

347 

0.006 

0.290 

0.006 

0.289 

0.208 

P79b 

Procubitermes 

388 

383 

1.97 

292 

0.0 

0.451 

0.0 

0.447 

0.203 

T29 

Thoracotermes 

98 

96 

1.96 

90 

0.033 

0.551 

0.033 

0.515 

0.303 

T82 

Thoracotermes 

1073 

1470 

2.74 

1069 

0.018 

0.177 

0.018 

0.164 

0.159 






mean 

0.033 

0.232 

0.032 

0.227 

0.216 



0/71 <j)/7C 


FIG. 1. (Color online) Empirical network description in the 
case of nest C9 [Cubitermes sp., Central African Republic, 
see Fig. IB in j28j]). (a) Each tunnel or edge is represented 
as a vector in spherical coordinates, (b) tunnel lengths (r) 
distribution, (c) distribution of the vertical component ( 9 , 0 
points upwards and n points downwards), and (d) distribution 
of the horizontal component (</>). (c) and (d) are symmetric 
and periodic respectively because both i — c 2 and 1?2 — c 1 
are used. 


ate a new node (and give it an increasing unique identi¬ 
fication number i > 1) or connect to an existing node if 
it is close to the end of the new edge, (v) prune edges of 
lesser importance. 

(i) Initial and boundary condition: From the em¬ 
pirical nest data we can compute the x-y-z inter¬ 
vals [x m im X rnax \i [//mm ? and [z rfl7n . Z m ax] for 

each nest. We impose initial and boundary condi¬ 
tions based on these intervals. We then choose an 
initial node (xo,yo, zo) randomly in the intervals 
Xo £ [0.75 Xmin + 0.25x maa ,, 0.25 Xmin 0 .75x maa ,j, 


2/o £ [0-75y m i n +0.25y max , 0.25y m j n -|-0.75// ma;c ] and 
z = Zmi n . During nest growth we assume boundary 
conditions which have the shape of an ellipse in the 
x — y plane. The formula of the ellipse is given by 


(^r) 


y-Y c 

B 


2 


= 1 


(1) 


with A c — ( X m i n T X m ax) /2, Y c — ( 2/mzn “b 2 Imax) 

A = ( Xmax Xmin) /2 and B = (l Imax //mm) /2. 
The z coordinate cannot be below z m i n , but there 
is no upper limit for z. The biological rationale 
behind these choices is that termite nest volumes 
are proportional to colony size [3ll . colony size is the 
number of individuals living in a single colony] and 
that nests are not enlarged during colony growth 
but rather rebuilt from scratch - colony size can 
therefore provide a template for the surface of the 
initial nest site construction and for the final nest 
height. 

(ii) Determining the node of the next edge’s ori¬ 
gin: At each time step we choose a node randomly 
from all existing ones. However, in order to fos¬ 
ter nest expansion (growth from peripheral nodes) 
we give preference to the latest added nodes by 
choosing node i with probability exp(—A(n — *)). 
Here n is the total number of nodes for the given 
nest. The characteristic time constant A is specific 
to each nest and will be chosen by calibration (see 
below). The biological idea is that termites mark 
construction material with some volatile chemical 
[32] and prefer recently built elements to continue 
construction. 

(iii) Direction and length of the next edge: The 

length (r) and the directions (0, (f >) of the next 
edge are chosen randomly from the normal or uni¬ 
form distributions determined above. Edge length 
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r is restricted to r > r m i n where r m i n is the ob¬ 
served minimum edge length (see also the next rule, 
shorter edges would lead to an edge connection to 
the original node). The vertical component 9 is re¬ 
stricted to the interval (0,7r). Furthermore, since 
nest construction goes upward and since edges in 
the original data are not oriented we avoid excessive 
downward construction by replacing a 8 > 7r/2 + ag 
by 9 = ir — 9 (recall that 9 = n points downwards). 
Finally, <f> is drawn randomly from a uniform dis¬ 
tribution in (0, 2n). 

(iv) Edge construction: This new edge is only added 
if it does not quit the ellipsoid boundary condi¬ 
tion (otherwise a new edge (r, 9, <j>) is drawn until 
it remains within these boundaries). The biologi¬ 
cal rationale is that termites sense gravity (33] and 
do not extend construction over empty space (note 
that the construction of the characteristic “hats” 
in Cubitermes nests is not included in our model). 
A new node is created if there is no existing node 
within the distance r m j n of the endpoint of the new 
edge. If there already exists a node within this dis¬ 
tance the new edge is connected to this node. 

(v) Pruning of edges: At each time step i we com¬ 
pute the edge-betweenness of each edge 0, that 
is the number of shortest paths between pairs of 
nodes that pass through this edge, and remove 
the one with the smallest value with probability 
£. Such pruning is a common feature of most ob¬ 
served transportation networks (reviewed in Qu¬ 
it has been directly observed in termite’s under¬ 
ground tunneling networks 0 m , and there is an 
indirect evidence that it happens in Cubitermes 
nest growth [28]. 

Repeat steps (ii-v) until the model has the same num¬ 
ber of nodes as the original nest. The biological ratio¬ 
nale underlying this stop criterion is again the observa¬ 
tion that mound size is proportional to colony size (31] 
and mound size is better approximated by the number of 
chambers (volume) than the number of edges (network 
length). 

Since A is related to nest height h and £ is related 
to the total number of tunnels L, we determine these 
two parameters one by one, first estimating A (with an 
arbitrarily fixed £) by minimizing the error function 

6l = | hdata h m odel\/hdata (2) 

with a standard bisection method, then estimating £ (fix¬ 
ing the already estimated A) by minimizing the error 
function 


to test whether the estimated parameters are specific for 
a genus, reporting the resulting F statistic with the as¬ 
sociated degrees of freedom and p-values. 


D. Alternative models: Random geometric graph 
(RGG) model and random pruning (RP) model 


We consider two alternative models to clarify the per¬ 
formance of the BBP model and the roles of local rules 
in the BBP model: (i) random geometric graph model 
and (ii) random pruning model. A random geometric 
graph (RGG) is the mathematically simplest spatial net¬ 
work [2j]. The RGG model for a given nest is embed¬ 
ded in the cylindrical space having the same interval 
of \Xmim %max\i [Urnim Umax]? and [C m7 ; n , Z rn rf , 7; ] obtained 
from the nest with the ellipsoid boundary condition given 
by Eq. (1). The nodes of the RGG model are uniformly 
distributed in this cylinder. Two nodes i and j are con¬ 
nected if the below condition is satisfied: 



Here X !$, 1), and Z, indicate x, y, and z coordi¬ 
nates of node i, respectively and A = ( x max — x m i n )/ 2 , 
R — (Umax brain ) /2 , and C — (^maa: Zmin )/2. We de¬ 
termine the value of Rrgg such that it gives us the same 
number of links in the RGG network as in the original 
nest network. 

The random pruning (RP) model is the same as the 
original model in every points (i.e., from (i) to (iv) in the 
previous subsection) except the priming process (i.e, (v) 
in the previous subsection) where we remove a randomly 
selected link with probability £ rather than the lowest 
edge-betweenness link. 


E. Model validation 

We will use five emergent properties to compare be¬ 
tween the original and the simulated networks: (i) node 
degree distribution, (ii) size of the largest connected com¬ 
ponent, (iii) average topological path length in the largest 
connected component, (iv) backbone link ratio (fraction 
of edges whose removal leads to a disconnection of the 
largest connected component), and (v) local graph re¬ 
dundancy (as defined in [36], it complements backbone 
link ratio by computing the mean of the inverse of the 
topological path length to connect two adjacent nodes 
once the direct link has been blocked: low values indi¬ 
cate long detours). Predicted properties will be based on 
1000 simulated networks. 


^2 — | Ldata Cmodel\/Ldata- ( 3 ) 

For each value of (A, £) we simulate 500 nests in order to 
compute mean h mo dei and L mo d e i • Overall error is de¬ 
fined as e = ei + 62 - We will use an analysis of variance 


F. Sensitivity analysis 

We assess the influence of A, £, and network size (num¬ 
ber of nodes) on four network properties: (i) size of the 
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largest connected component, (ii) average topological dis¬ 
tance between any two nodes of the largest connected 
network component, (iii) backbone link ratio, and (iv) 
local graph redundancy. For an efficient transport in¬ 
side the nest and easy defense against predators (such 
as ants) termite nests should show low average distance 
(fast transport), a high backbone link ratio (tree-like 
structures, tunnel blocking by a soldier efficiently isolates 
a part of the nest), and a low local graph redundancy 
(blocking of a tunnel forces attacking ants to take long 
detours). The values chosen for A and £ cover the range 
of the estimated values (from 0 to 0.4 in steps of 0.01), 
while the number of nodes cover the original nest sizes 
(200, 400, 600 and 1000 nodes). 


III. RESULTS 
A. Model validation 

Figure [2] shows the degree distributions of the real net¬ 
works and the simulated networks generated by the BBP 
model, RGG model, and RP model for each nest. The 
real networks and the simulated networks show qualita¬ 
tively similar behaviors, with peaks around k = 1,2 and 
exponential type decay of the right tail. However, all of 
the nests have a peak at k = 1 while the BBP model 
creates a peak at k = 2 for seven nests: the BBP model 
seems to create less dead ends than there exist in the real 
networks (see also Fig. Nest P67 shows the worst fit 
with a larger tail in the model networks: note that this 
nest looks like two nest parts that have fused together 
during growth (see SM Fig. 6), a process not considered 
in the simulated networks. The degree distributions gen¬ 
erated from the RGG model and the RP model show sim¬ 
ilar patterns as the distributions from the BBP model. 

The largest connected component (LCC) of the net¬ 
work is fundamental for internal transportation since 
communication is not possible between disconnected 
components. Figure @|a) represents the sizes of the 
largest connected components in real networks and sim¬ 
ulated networks. The BBP model generated networks 
having comparable sizes of the largest connected com¬ 
ponents with the ones in real nests for most cases while 
the RGG model failed to generate the largest connected 
components of proper size in most cases. However, the 
performance of the RP model is comparable with the 
BBP model. It is notable that the discrepancy between 
real networks and networks generated by RP models is 
larger when the average degree is low. The blind cutting 
of random pruning can in this case increasingly affect 
important links and thus reduce LCC, while with higher 
average degrees many ’back-up’ links exist that help pre¬ 
serve the LCC. We can observe this pattern also in the 
sensitivity analysis of the RP model (see Fig. 13 in 0): 
with increasing probability to prune an edge (£) LCC 
quickly degrades, while it is better preserved in the BBP 
model (see Fig.[5]below). This indicates that betweenness 


based pruning is a better strategy than random pruning 
to keep the largest connected component. 

Figure SKb) compares the average topological distance 
between any two nodes in the largest connected compo¬ 
nent, Fig. |4|c) the backbone link ratios, and Fig. |4](d) 
the local graph redundancy. The RGG statistics are 
not shown because computing them makes little sense 
when the LCC is too small. Overall both the BBP and 
RP model successfully reproduced the average distances, 
backbone link ratios, and local graph redundancies of 
the real networks. However, in the case of average dis¬ 
tance neither BBP nor RP catch the variation between 
nests: Kendall’s correlation coefficient r with the real 
nests ( tbbp = 0.30 and trp = 0.363) is not significantly 
different from 0. In the case of LCC (tbbp = 0.939 and 
prp = 0.818), backbone link ratio (r = 0.606 for both 
models) or local graph redundancy ( tbbp = 0.636 and 
t rp = 0.606) t is significantly different from 0. Given the 
plotted standard deviations notable differences can only 
be detected for nests C9, P67, P78, and T29 concerning 
average distance, nests CIO and P67 concerning back¬ 
bone link ratio, and CIO and P67 concerning local graph 
redundancy. Again, nest P67 stands out in this compar¬ 
ison. See [29| for the full distributions of distances as in 
Fig. 0 

To compare the performances of the BBP model and 
the RP model quantitatively we define a Z-score for each 
metric X such that: 

ry [ X rea l — (X) Model 

Z{X) = - (5) 

& Model 

where (X) Model is the average of X for the given model 
and a Model is the standard deviation. We show the 
Z-scores for each metric and each nest in Tables I-IV 
in 01- Interestingly we found that the BBP model has 
a tendency to perform better than the RP model for Cu- 
bitermes and Thoracotermes nests while for the nests of 
Procubitermes the RP model performs better than the 
BBP model. 


B. Calibrated parameters 

For each nest (except P67 and P79b) we found pa¬ 
rameter sets (A bbp^bbp) with e < 0.025 for the BBP 
models and ( Xrp , £rp) with e < 0.027 for the RP model. 
Nest P67 with cbbp = 0.285 and epp = 0.283 is again 
an exception. In the case of Nest P79b, the BBP model 
( £bbp = 0.125) and the RP model (crp = 0.129) pro¬ 
vide us the taller nests than real ones since the errors e 
are mainly from the height of the models (i.e., from ei in 
Eq. [2]) . Table [T] summarizes the principal nest properties 
and the estimated parameter sets. 

Model parameter £ does not depend on taxonomy 
(analysis of variance at the genus level, F = 0.78, df = 
(2,9 ),p = 0.49 for BBP and F = 0.60, df = (2,9 ),p = 
0.57 for RP), but A depends on Genus (F = 6.61, df = 
(2,9 ),p = 0.017 for BBP and F = 6.16, df = (2,9 ),p = 



6 



0 2 4 6 8 10 12 14 
k 


0 2 4 6 8 10 12 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 

k k k 



Q_ 



CL 




0 2 4 6 8 1012141618 
k 



0 2 4 6 8 10 12 14 


02468 10121416 
k 



FIG. 2. (Color online) Degree distribution of the real networks (Data: red line) and the networks generated by the BBP model 
(BBP: black line), by the RGG model (RGG: blue line), and RP model (RP: magenta line). P(k) are relative frequencies 
(summing to 1). 


0.021 for RP; Tukey post-hoc: Cubitermes has signifi¬ 
cantly higher A than Procubitermes and Thoracotermes). 

C. Sensitivity to parameters (A, £) 

Figure [5] shows the sensitivity of network properties 
to the two free parameters (A bbp^bbp)- We can ob¬ 
serve that for the given parameter space, the size of the 
largest connected component is preserved when £ < 0.3. 
We further see that average distance only depends on 
the decay rate A, with low distances for low values of 
A. Note that both sensitivities are quite different in the 
BP model (Fig. 13 in M)- the network quickly breaks 
down with increasing £, especially for large A, leading to 
small LCC’s and low average distances in these LCC’s. 
The nests have actually all values of A < 0.1, indicat¬ 
ing a moderate preference to continue construction from 
the most recent nodes. Backbone link ratio and local 
graph redundancy have a more complex dependence on 
(A,£) (isolines seem to be linked to the product A£), but 
termites could obtain a high backbone link ratio and a 
low local graph redundancy by increasing both A and £. 


All these qualitative observations are independent of nest 
size (number of nodes). The fact that the analyzed nests 
have nevertheless A < 0.1 indicates that low average dis¬ 
tances have more importance than increasing backbone 
link ratio or decreasing local graph redundancy. Only in¬ 
creasing £ would further optimize the last three criteria, 
but also produce more disconnected chambers (decrease 
LCC) which incurs a cost to the colony (construct living 
space that cannot be used). 

IV. DISCUSSION 

In this paper we develop a simple network growth 
model (the BBP model) to test how the nest architecture 
of termites in the Termitinae sub-family emerges from 
self-organization based on local rules only. These nests 
consist of spherical chambers connected by tunnels, an 
architecture that can be represented as a network. The 
BBP model uses empirical nest descriptions (edge length 
and orientation) and two free parameters that control pe¬ 
ripheral growth and pruning of existing edges. The free 
parameters are calibrated by fitting the model to the ob- 
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FIG. 3. (Color online) Analyzed Thoracotermes nest T29 
with its network representation (middle) and a simulated net¬ 
work (right) by the BBP model. The pink node (lower right) 
indicates the initial node of the simulated network. Note that 
while there are isolated chambers in the original network, iso¬ 
lated chambers are more prevalent in the simulated nests. See 
the Discussion for further comments. 


served nest height and the number of edges. 

The BBP model correctly reproduces several emer¬ 
gent properties: a) node degree distribution (Fig. [2]), b) 
the size of the largest connected component (Fig. 0(a)), 
c) the average distance between any two nodes in the 
largest connected component (Fig. 0(b)), d) the back¬ 
bone link ratio (Fig. 0(c)), and e) the local graph redun¬ 
dancy (Fig. 0(d)). The properties (c) to (e) are of eco¬ 
logical relevance: the populations move around in these 
nests to place brood or larvae in chambers with optimal 
climatic conditions or to store and retrieve food (dead or¬ 
ganic matter), short distances are therefore useful. It has 
already been shown for Cubitermes nests that their av¬ 
erage distances are shorter than what could be obtained 
by randomly connecting the existing nodes (36| . and our 
model correctly reproduces these lengths for most an¬ 
alyzed nests. The nest must also protect the colony 
against ant predators: this is done by the soldiers who 
can block a tunnel with their head capsule. A treelike 
network structure (or high backbone link ratio combined 
with local graph redundancy) helps with this strategy, 
and our model also correctly predicts these properties. 

We compare our BBP model to two alternatives, the 


RP model that is identical to BBP but with pruning ap¬ 
plied to randomly chosen edges instead of the smallest 
betweenness centrality edge as in BBP, and to the Ran¬ 
dom Geometric Graph (RGG) model with the same space 
constraints as the two previous models. The latter can 
be ruled out because it cannot reproduce the observed 
largest connected component (LCC). The RP model gen¬ 
erally performs as well as the BBP model (in the case of 
Procubitermes even better according to the Z-scores) for 
nests with average degree >2.0. However, in nests with 
lower average degree (C10,Cll,P79b,T29) RP tends to 
give fragmented networks. Furthermore, the sensitivity 
analysis shows that LCC size is conserved in the BBP 
model for a much larger range of pruning probabilities 
than in the RP model. Overall, pruning based on be¬ 
tweenness centrality (BBP) reproduces networks similar 
to the real networks more robustly than the RP model, 
we therefore conclude that betweenness centrality based 
pruning might be an important mechanism in termite 
nest construction. 

While our model is intended to specifically reproduce 
the general properties and appearance of termite nests of 
the Termitinae sub-family, we can speculate that the key 
ingredients that regulate network morphogenesis in our 
model (peripheral growth and pruning) are likely to be 
shared also by a large number of biological transporta¬ 
tion networks. For instance, pruning phenomena are ob¬ 
served in the maturation of neural networks (e.g. thr oug h 
programmed cell death [37j and synaptic pruning [381]). 
Mycelial networks formed by fungi similarly undergo a 
maturation process that involves peripheral growth of the 
hyphal tips in response to local changes of turgor pres¬ 
sure [39| and regression of filaments from nutrient de¬ 
pleted regions [40(. Growth of peripheral filaments that 
are subsequently pruned is also central to the formation 
of the network of cytoplasmic tubes that constitute the 
body of the plasmodium of Physarum polycephalum fill . 
It seems plausible that self-organized network construc¬ 
tion and optimization requires these mechanisms to op¬ 
erate on an initially highly connected network that is 
subsequently pruned. This could respond for instance to 
the fact that self-organized mechanisms can only evalu¬ 
ate transportation performance through the transporta¬ 
tion itself, and the fine-tuning and optimization of the 
network requires to operate on already formed connec¬ 
tions. 

The sensitivity analysis (Fig. 0]) indicates that our two 
free parameters A (for peripheral growth) and £ (for prun¬ 
ing) strongly influence both average distance and back¬ 
bone link ratio (as well as local graph redundancy, an¬ 
other measure of how efficiently the nest can be defended, 
see jH|). We calibrated them by concentrating on a mini¬ 
mal number of features, nest height and total number of 
tunnels (i.e., links), in order to use the other nest fea¬ 
tures as emergent properties for model validation. The 
good qualitative agreement between these emergent fea¬ 
tures and the original nest features indicates that the 
estimated parameter values would not change much if 
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FIG. 4. (Color online) Comparison between the real networks (Data: red bar) and the networks generated by the BBP model 
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Local graph redundancy. The error bars represent standard deviations computed from 1000 model generated networks. 


calibration had been based on more features. The pe¬ 
ripheral growth parameter A bbp has mean value 0.033, 
meaning that the probability for a given node to be cho¬ 
sen for the next edge is divided by half every 30 newly 
added chambers. Typical Cubitermes nests grow in 1-3 
months [30j, meaning that the half life is of the order of 
some days to a week, a plausible duration for chemical 
marking in termites. In larger nests our model will also 
predict that growth becomes spotted on the nest surface, 
in agreement with field observations by one of the authors 
(CJ, unpublished). 

Our model does not explain how termites decide to dig 
a new tunnel or how they choose an orientation in space 
- it is an empirical model at an intermediate scale. As 
such it resembles the 2D ant tunneling network models 
suggested by @,d, ISL or 2D termite tunneling mod¬ 
els as suggested by @-Q. However, though 3D network 
data of social insect nests become increasingly available 
[43]-|45j], our model seems to be the first to predict the 


nest’s 3D network architecture. Further work should in¬ 
vestigate how nests grow in time [45 ] and, on the other 
end, how termites decide where and when to construct. 

Another feature not predicted by our model are the 
characteristic “hats” atop the Cubitermes nests (see SM 
[29j| ). It is not known how termites decide to start con¬ 
structing laterally, and we could not identify statistical 
properties specific to the height where the hat is con¬ 
structed. The behavioral algorithm underlying this spe¬ 
cific “hat” feature is therefore an open question. 

There were also some misfits in this work. For example, 
nest P67 is badly explained. This nest actually consists 
of two columns that have grown together (see SM [29l|). 
This bad fit might be corrected by letting the model grow 
from two randomly chosen initial nodes. Also, the above 
ground nest is only part of a termite colony, it connects 
to an extensive underground tunneling network through 
which termites access food. The absence of this under¬ 
ground network in our data leads to an underestimation 
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of edge-betweenness of the lower edges, thus explaining 
the often observed priming of these edges (Fig. [3]) that 
leads to an illogical disconnection between the above and 
below ground nest parts. This could only be corrected if 
one collects, in addition to the above ground par t, a cast 
of the corresponding underground network [44J | . Note 
also that underground tunnels are often built at a con¬ 
stant distance below ground, thus effectively leading to a 
2D below ground network [2|| |46[ : an extended dataset 
could therefore explore how 3D networks connect to 2D 
networks and how this alters network properties. A fur¬ 
ther discrepancy is observed in the node degree distri¬ 
butions (Fig. [2]): all nests have a peak at degree 1, but 
for half the nests the model predicts a peak at degree 2. 
This could probably be corrected by including this cri¬ 
terion when fitting A, or by treating the fusion distance 
R (model part iv) as a free parameter. In the interest 
of keeping the model simple we refrained from such ex¬ 
tended fitting procedures. 

The sensitivity analysis (Fig. [5]) shows that average 
distance only depends on the peripheral growth param¬ 
eter A, not on the pruning parameter £. Backbone link 
ratio and local graph redundancy, on the other side, de¬ 
pend on both parameters (higher A can be compensated 
by lower £ or the other way around). Both observations 
are true for all analyzed nest sizes. An efficient termite 
nest should have small average distances, high backbone 
link ratios and low local graph redundancy. The detected 
nest positions suggest that short distances are more im¬ 


portant for the colony than the other two criteria (that 
can be linked to nest defense). However, further infor¬ 
mation is required on what termites actually do in their 
nests and on their vulnerability to predation before fur¬ 
ther speculating about such issues. 

In sum, we found a parsimonious empirical network 
growth model based on self-organized principles that 
successfully predicts the nest architecture of Termitinae 
nests. Peripheral growth (i.e. some volatile chemical 
marking of new chambers/nodes/edges) and pruning of 
less important edges are important ingredients in this 
model. 
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